function modelResults()

  clc
  addpath ./utils/ -end

  N=10;
  %step=1e-4;
  step=1e-3;
  errTol=step/2;
  vLambda=[0.001 0.002 0.004 0.006 0.008  0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.15 0.2 0.3 0.4 0.6 0.8 0.9 1.0];

  %% Fig 5, 6 and 7
  %typeStr={'$CW^{n1}=1$', '$CW^{n2}=2$'};         
  typeStr={'CW^{n1}=1', 'CW^{n2}=2'};         
  CW    =[2,    2];
  maxBOs=[4,    4];
  aMinBE=[3     3];
  aMaxBE=[5,    5];
  vM    =[12,   0];  %Simple way to simulate nodes of just one type ...

  [resCMatrix, indRV]=solveModel_12_2(CW, maxBOs, aMinBE, aMaxBE, vM, N, vLambda, step, errTol);
  titleStr=['CW=[', num2str(CW), '], maxBOs=[', num2str(maxBOs), '], aMinBE=[', num2str(aMinBE), '], aMaxBE=[', num2str(aMaxBE), ']'];
  plotResults(titleStr, typeStr, resCMatrix, indRV);

  %% Fig 8 and 9
  %typeStr={'$CW^{n1}=1$', '$CW^{n2}=2$'};         
  typeStr={'CW^{n1}=1', 'CW^{n2}=2'};         
  CW    =[1,    2];
  maxBOs=[4,    4];
  aMinBE=[3     3];
  aMaxBE=[5,    5];
  vM    =[12/2, 12/2];

  [resCMatrix, indRV]=solveModel_12_2(CW, maxBOs, aMinBE, aMaxBE, vM, N, vLambda, step, errTol);
  titleStr=['CW=[', num2str(CW), '], maxBOs=[', num2str(maxBOs), '], aMinBE=[', num2str(aMinBE), '], aMaxBE=[', num2str(aMaxBE), ']'];
  plotResults(titleStr, typeStr, resCMatrix, indRV);

  %% Fig 10 and 11
  %typeStr={'$M_{bo}^{n1}=1$', '$M_{bo}^{n2}=5$'};                                                       
  typeStr={'M_{bo}^{n1}=1', 'M_{bo}^{n2}=5'};                                                       
  CW    =[2,    2];
  maxBOs=[1,    5];
  aMinBE=[3     3];
  aMaxBE=[5,    5];
  vM    =[12/2, 12/2];

  [resCMatrix, indRV]=solveModel_12_2(CW, maxBOs, aMinBE, aMaxBE, vM, N, vLambda, step, errTol);
  titleStr=['CW=[', num2str(CW), '], maxBOs=[', num2str(maxBOs), '], aMinBE=[', num2str(aMinBE), '], aMaxBE=[', num2str(aMaxBE), ']'];
  plotResults(titleStr, typeStr, resCMatrix, indRV);

  %% Fig 12 and 13
  %typeStr={'$\mathbf{p}^{n1}_{b} [0,5]$', '$\mathbf{p}^{n2}_{b} [3,5]$'};                               
  typeStr={'p^{n1}_{b} [0,5]', 'p^{n2}_{b} [3,5]'};                               
  CW    =[2,    2];
  maxBOs=[5,    5];
  aMinBE=[0     3];
  aMaxBE=[5,    5];
  vM    =[12/2, 12/2];
  
  [resCMatrix, indRV]=solveModel_12_2(CW, maxBOs, aMinBE, aMaxBE, vM, N, vLambda, step, errTol);
  titleStr=['CW=[', num2str(CW), '], maxBOs=[', num2str(maxBOs), '], aMinBE=[', num2str(aMinBE), '], aMaxBE=[', num2str(aMaxBE), ']'];
  plotResults(titleStr, typeStr, resCMatrix, indRV);

  %% Fig 14 y 15
  typeStr={'$CW^{n1}=1, \mathbf{p}^{n1}_{b} [0,5]$', '$CW^{n2}=2, \mathbf{p}^{n2}_{b} [3,5]$'};         
  CW    =[1,    2];
  maxBOs=[5,    5];
  aMinBE=[0     3];
  aMaxBE=[5,    5];
  vM    =[12/2, 12/2];

  [resCMatrix, indRV]=solveModel_12_2(CW, maxBOs, aMinBE, aMaxBE, vM, N, vLambda, step, errTol);
  titleStr=['CW=[', num2str(CW), '], maxBOs=[', num2str(maxBOs), '], aMinBE=[', num2str(aMinBE), '], aMaxBE=[', num2str(aMaxBE), ']'];
  plotResults(titleStr, typeStr, resCMatrix, indRV);

end
  
function plotResults(titleStr, typeStr, resCMatrix, indRV)
  
  lambda=cell2mat(resCMatrix(:,indRV.lambda));
  
  %% Throughtput
  mS=cell2mat(resCMatrix(:,indRV.S));
  legStr={};  
  figure;
  loglog(lambda, mS(:,1),  'k-');   legStr{end+1}='Aggregate';  hold on;
  loglog(lambda, mS(:,2),  'b-');   legStr{end+1}=typeStr{1};
  loglog(lambda, mS(:,3),  'r-');   legStr{end+1}=typeStr{2};
  title(['Throughtput for ', titleStr]);
  %legend(legStr, 'Location', 'NorthWest','interpreter', 'latex')
  legend(legStr, 'Location', 'NorthWest')
  xlabel('\lambda [Packets per Packet Duration]'); ylabel('Throughtput');
  axis([1e-3 1e0 1e-2 1e0]);
  loglog([0.05, 0.05],  [0.01 1], 'k');
  set(gcf, 'PaperPositionMode', 'auto');

  %% Latency
  mLatency=cell2mat(resCMatrix(:,indRV.latency));
  legStr={};  
  figure;
  semilogx(lambda, mLatency(:,1),                            'Color', 'r');                      legStr{end+1}=['Lat. ',     typeStr{1}]; hold on;
  semilogx(lambda, mLatency(:,2),                            'Color', 'r', 'LineStyle', '--');   legStr{end+1}=['Lat. ',     typeStr{2}];
  semilogx([0.05, 0.05],  [0.01 max(max(mLatency))], 'k');
  title(['Latency for ', titleStr]);
  %legend(legStr, 'Location', 'NorthWest','interpreter', 'latex')
  legend(legStr, 'Location', 'NorthWest')
  xlabel('\lambda [Packets per Packet Duration]'); ylabel('Latency [BackOff Slots]'); ylabel('Latency [BackOff Slots]');
  
  %% Probability of Delivery
  pReject=cell2mat(resCMatrix(:,indRV.pReject));
  pSend=cell2mat(resCMatrix(:,indRV.pSend));
  mPDR=cell2mat(resCMatrix(:,indRV.PDR));  
  pDeliver=[(1-pReject(:,1)).*pSend(:,1).*mPDR(:,2) , (1-pReject(:,2)).*pSend(:,2).*mPDR(:,3)];
  legStr={};  
  figure
  semilogx(lambda, pDeliver(:,1),  'Color', 'b');                      legStr{end+1}=['P. Del. ', typeStr{1}]; hold on;
  semilogx(lambda, pDeliver(:,2),  'Color', 'b', 'LineStyle', '--');   legStr{end+1}=['P. Del. ', typeStr{2}];
  semilogx([0.05, 0.05],  [0.01 max(max(pDeliver))], 'k');
  title(['P. Deliver for ', titleStr]);
  %legend(legStr, 'Location', 'SouthWest','interpreter', 'latex')
  legend(legStr, 'Location', 'SouthWest')
  xlabel('\lambda [Packets per Packet Duration]'); ylabel('Latency [BackOff Slots]'); ylabel('Probability');

  %% Power Consumption
  vWRx=cell2mat(resCMatrix(:,indRV.WRx));
  vWTx=cell2mat(resCMatrix(:,indRV.WTx));
  vWIdle=cell2mat(resCMatrix(:,indRV.WIdle));
  vWTotal=vWRx+vWTx+vWIdle;

  legStr={};  
  figure;
  loglog(lambda, vWTotal(:,1),  'b');   legStr{end+1}=['Total Power ', typeStr{1}]; hold on;
  loglog(lambda, vWTotal(:,2),  'r');   legStr{end+1}=['Total Power ', typeStr{2}];
  loglog(lambda, vWRx(:,1),     'b-.'); legStr{end+1}=['Rx Power ',    typeStr{1}];
  loglog(lambda, vWRx(:,2),     'r-.'); legStr{end+1}=['Rx Power ',    typeStr{2}];
  loglog(lambda, vWTx(:,1),     'b-*'); legStr{end+1}=['Tx Power ',    typeStr{1}];
  loglog(lambda, vWTx(:,2),     'r-*'); legStr{end+1}=['Tx Power ',    typeStr{2}];
  loglog(lambda, vWIdle(:,1),   'b-o'); legStr{end+1}=['Idle Power ',  typeStr{1}];
  loglog(lambda, vWIdle(:,2),   'r-o'); legStr{end+1}=['Idle Power ',  typeStr{2}];
  %legend(legStr, 'Location', 'NorthWest','interpreter', 'latex')
  legend(legStr, 'Location', 'NorthWest')
  title(['Power Consumpion for ', titleStr]);
  xlabel('\lambda [Packets per Packet Duration]'); ylabel('Power [mW]');
  axis([1e-3 1e0 1e-2 1e2]);

end
